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1. Introduction 

The electronic properties of graphene - a two-dimensional monolayer of graphite 
made of carbon atoms only - have recently attracted a lot of interest [U El E] . For 
energies close to the charge neutrality point, noninteracting quasi-particles in graphene 
(henceforth called "electrons") are well described by the Dirac- Weyl Hamiltonian of 
massless relativistic fermions. This suggests an easily accessible condensed-matter 
realization of relativistic quantum mechanics. Recent interest has turned to Coulomb 
interaction effects, in particular to the case when N electrons are confined to a 
finite region (a so-called "quantum dot" ) in the graphene layer. Using electrostatic 
confinement potentials, such Coulomb-correlated artificial atoms have been studied 
in detail for the two-dimensional electron gas in semiconductors @], where the non- 
relativistic Schrodinger equation describes the single-particle sector. In graphene, 
however, the Klein tunnelling phenomenon [2] of relativistic Dirac fermions renders the 
standard electrostatic confinement method experimentally difficult or even impossible 
to apply, and usually only resonances can be expected [2 HI |9] . As an alternative, 
confinement by inhomogeneous orbital magnetic fields has been suggested [10] . and 
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the electronic structure of two interacting electrons in such a magnetic quantum dot 
in graphene was recently studied using exact diagonalization [TT] . Experimentally, 
up to now, lithographically fabricated quantum dots were mostly studied |12) . where 
the boundary is rather disordered, the confinement potential cannot be tuned, and 
a detailed comparison of experimental data to theory is difficult. On the other 
hand, inhomogeneous magnetic fields have been experimentally generated and studied 
in semiconductor devices by using suitable lithographically deposited ferromagnetic 
layers [13], and the generalization to graphene should pose no fundamental obstacle. 
Concrete experimentally relevant profiles for quantum dot confinement by such fields 
were theoretically studied also in Ref. [14] . We note that an artificial vector potential 
giving rise to the same mathematical model can also be generated by applying 
mechanical forces, producing appropriate strain in the sample [2] [3]. 

In this paper, given the widespread interest in understanding and usefully 
employing the electronic structure of graphene quantum dots, we address the definition 
and the stability of the relativistic interacting iV-particle system in such a magnetic 
graphene dot, primarily from a mathematically oriented perspective. To that end, we 
analyze in detail the Hartree-Fock functional and show that under certain conditions, 
a minimizer exists. The maximum number N c of particles is computed and shown 
to depend on the interaction strength. When N < N c we solve the Hartrcc-Fock 
equations numerically. For N > N c , no minimizer with particle number N exists. 
The excess particles are not localized and appear, numerically, occupying bulk Landau 
states centered far away from the dot region to lower the repulsive interaction energy. 

It is well known that for relativistic iV-particle problems, the presence of 
interactions implies that one should use a projection scheme (T5j HH [17], whose 
integrity for graphene dots - despite the fact that bulk graphene corresponds to a 
gapless model - has been shown in Ref. [11] . At first sight, an inclusion of the entire 
Dirac sea along the more fundamental lines of Hainzl et al. (see [18] and the references 
therein) might seem desirable. However, this would complicate matters unduly, since 
- after all - the description of graphene by a two-dimensional Dirac equation emerges 
from a non-relativistic band structure calculation. Moreover, in the presence of a 
finite gap separating occupied and empty states - which is the case below - it is not 
only reasonable to freeze the Dirac sea as given by the external field (Furry picture), 
but also to assume that no electron-positron pairs are created. It should be remarked 
that this strategy is not only supported by Ref. pj] but is a standard procedure in 
quantum chemistry [19] . Because of this we will restrict ourselves in this work to the 
no-pair model, more precisely to the no-pair model in the Furry picture with given 
electron number N . 

In this model, the electronic states of graphene in the presence of a magnetic field 
V x A are the unit vectors in X(o,oo)(«<x • (— iftV + ~ A))[L 2 (1R 2 : C 2 )], where A is the 
magnetic vector potential, v is the Fermi velocity in graphene, cr = (01,02) are the 
first two Pauli matrices, — e < the charge of the electron, and h the Planck constant. 
We introduce convenient units by scaling momentum p t— > p/(vh) and coordinates 
x 1 y vhx, which generates a unitary transform. N Coulomb-interacting electrons in 
graphene are then formally described by the Hamiltonian 

N 

k • (p« + A ( x «)) - ^( x »)] + 1 — z — 1 w 

n=l l< m <n<N |Xm Xn| 

projected onto the antisymmetric tensor product of the above space. Here A(x) := 
(ev I c)A(vhx.) and we have added an external electrostatic potential ip. The interaction 
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strength is encoded in the dimensionless fine structure constant a = e 2 /(vnh). 
Physical values in graphene are < a < 2, with the precise value depending on 
the dielectric constant k of the environment (e.g., the substrate on which graphene 
is deposited). The upper limit for a is approached only for suspended samples. It 
is well known that quasi-particles in graphene have additional valley ( U K point") 
and electronic spin degrees of freedom. For field configurations that are smooth on 
the scale of grapheme's lattice constant ao = 0.246 nm, however, no valley mixing is 
expected, and our simpler description with just one K point and full spin polarization 
is sufficient This is the situation considered in our paper. In any case, the 

generalizations necessary to go beyond the single-spin and single-valley model Jl]) are 
conceptually unproblematic. A similar approach is also expected to apply for magnetic 
dots in bilayer graphene, where the linear dispersion relation of monolayer graphene 
is modified [2J. 

The outline of the remainder of this paper is as follows. In Sec. [2J we specify and 
explain the mathematical model. The Hartrec-Fock functional is considered in Sec. (3J 
and we prove the existence of a minimizer. For the benefit of the mathematically 
oriented reader, we have included all proofs in detail. To make the abstract discussion 
in Scc.[3]concrete, we describe a specific magnetic field profile A(x) in ScclU modelling 
a circular magnetic dot of radius R around the origin. The confinement is here 
generated by a constant magnetic field B for r > R but zero field for r < R. The 
relativistic interacting A-particle problem for a magnetic dot in graphene is then 
studied for this field profile in Sec.|4]by numerically solving the Hartree-Fock equations 
discussed in Sec. [3] This calculation also yields quantitative predictions concerning 
the stability of the iV-electron system. We conclude with an outlook in Sec. El Some 
technical details have been relegated to two Appendices. 

2. Model 

To describe the physical situation where electrons are bound to a magnetic quantum 
dot, the magnetic vector potential A = Ao + Ai is expressed as the sum of a 
homogeneous magnetic field of strength B > perpendicular to the graphene plane, 
Ao(x) = (B/2)(— X2, Xi), and a perturbation Ai. For total external electromagnetic 
field A := (<p, A), using $ := cr ■ (p + Ao), we define 

D A :=p-<p:=a-(p + A)-ip = p Q + P (2) 

with the "perturbation" P := cr ■ Ai — (p. We define the one-electron Hilbert space 
with respect to A as 

Sj A := A+[L 2 (R 2 : C 2 )] , A+ := X( o,oc) (Da) ■ (3) 

Although more general cases can be treated, we assume for simplicity that P is a 
bounded operator which is sufficient for the application discussed and minimizes 
the amount of technical arguments needed. In the same spirit, we require relative 
compactness of the perturbation^ 

PeW and \P\i\$ Q -Li\-^ e&°°(Sj A ) ■ (4) 

Physically, this means that the perturbing electromagnetic potential decays at infinity 
and has only moderate singularities whose Fourier coefficients can be controlled by 

t For p S [l,oo), we set 6 p (^a) = {A & B(f) A )\ tr \A\ P < oo}, write e°°(#A) for the space of 
compact, and 'S(S^a) for the bounded operators on S)a- 



Multiparticle equations for graphene 



4 



the kinetic energy. At the same time, P is responsible for creating bound states that 
define the quantum dot. 

The energy of an iV-particlc state ip S 5(R 2Ar : C 2 ™) n -f^H is then given by 



Ety) := V, 



N 



n=l l<m<n<N 



if, . (5) 



Throughout this paper, fx is a positive constant smaller than the first positive 
eigenvalue (first bulk Landau level) of p . Of course, since \i is just a constant, the 
energy is shifted merely by —fitrj. However, this shift serves an important technical 
purpose: it will allow us to replace the minimization under the constraint trj < N 
instead of trj = N. If the minimizer 7 has trace N, then the problem with the N 
electron constraint has been solved. Note, that it may happen that the trace of the 
minimizer, tr7, stays always below some n c , even if fi is close to the Landau level. 
Then the smallest such n c - call it iV c - is the maximal number of electrons which can 
be captured by the dot, see Sec. [4] For N > N c , there are N — N c unbound electrons. 
The case that there is no minimizer with trace equal to N corresponds exactly to the 
situation, where those N — N c electrons float away to infinity and cannot be bound 
by the dot. 

The form ([3]) is obviously bounded from below and closable. This allows us to 
define the Hamiltonian Ba as its Fricdrichs extension. The Friedrichs extension is the 
canonical way to construct a self-adjoint Hamiltonian out of a symmetric operator 
which is bounded from below, see, e.g., Refs. [2Ql Satz 4.15] or [21] Theorem X.23]. 



3. The Hartree-Fock Functional 

3.1. The Hartree-Fock Functional 

We are now interested in the Hartree-Fock approximation for the ground state of the 
Hamiltonian Ba- The Hartree-Fock method is a standard tool to access interaction 
effects in atomic, molecular, or condensed- matter systems P33 The Hartree- 

Fock ground-state energy provides an upper bound for the true ground-state energy 
corresponding to Eq. (|5|), and it can be used to assess the stability of the iV-particle 
problem (23] . We mention that a different variational approach based on the so-called 
Muller functional can yield lower bounds for the ground-state energy [24j [25], and 
very useful results can be obtained by combining both methods. 

With the above choice for fi, let us denote by d the distance of /1 to the nearest 
spectral point of p . We start our analysis by defining the basic class of operators (7) 
that enter the Hartree-Fock functional. In physical terms, 7 is the density operator. 

Definition 1 We define the Banach space 

F:={je <B(S a )||| 7 ||f := \\\p - v\ 1/2 j\p - m| 1/2 ||i < 00, 7 = 7*} 

Here, as customary, ||a||i := tr y/a*a denotes the trace norm of the operator a. This 
definition of F is motivated[j] by the fact that the state represented by 7 should have 

§ Here S denotes the space of Schwartz functions, and f)^ is the TV-fold antisymmetric tensor 
product of the one-electron Hilbert space fiAi i- e -> the canonical N electron space. 
I For later use of the Banach-Alaoglu theorem, we also note that 

F, :={*6®(«A)|||» -Ml~ 4 «llto-Ml~* ee°°(f,A)} 
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finite kinetic energy, finite particle number, and real occupation numbers, the latter 
being the eigenvalues of 7. Note that S^-f)^) D F since 



i\\f = tr J\p Q -Ml 2 7l^ -Vhlpo - Ml 



>d»ta^|j» -/*|*|7||7||j» -/i|3 

= d* tr ^|7|||* -Hl7l > dix | 7 | ■ 

For a given element 8 G F, we denote its eigenvalues by A„ and its eigenspinors by £„. 
The associated integral kernel 8{x,y) if0 

^ y) : = A «^™( a; )^(y) • ( 6 ) 

n 

Associated with 8 is its one-particle density 

2 

WW :=££A„|£„(*)| 2 , 

s=l n 

its electric potential ip^ (x) := J dy ^(y)/|x — y|, and its exchange operator X^ s > in 
terms of the integral kernel y) = 5(x, y)|x - y| _1 . Then = (p {s) - is 

the mean field potential, and the Weyl operator associated to 8 is 

h { H :=D A -fi+ aWW . (7) 

The Coulomb scalar product is defined as 



D M :=U dx/ dy^^, (8) 
2 Jr2 J M 2 |x - y| 

and the exchange scalar product for 7,7'eF is 

Indicating quadratic forms of sesquilinear forms by brackets, e.g., X(j, 7) = X[y], we 
setQ[ 7 ] :=£>K]-X[ 7 ]. 

We now discuss inequalities which will allow us to define the Hartree-Fock energy 
functional. 

Lemma 1 We have 

|p + Ao| < |pj + B 1 ' 2 < \p Q - fx\ + p + B 1 ! 2 . 

Proof. We have 

(p + A ) 2 V) < (i>, ((P + A ) 2 + <x B + = W>, (f + B)xj>) . 



Using that the root is operator monotone and yp 2 , + B < I^J + B 1 / 2 , the first 
inequality follows. The second inequality is then clear. □ 

is a Banach space for which F is the dual space. The duality is given naturally as 

(7,5) :=tr(|j4 -/x|57| } 4 -M|5| J 4 - M |-|a|ji - /[ x|-3^ = tr( 7 <5) . 

We use the notation x = (x, s) for an element of G := K 2 X {1, 2} and dx for the product of the 
Lebesgue measure on R 2 with the counting measure in {1, 2}. 
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Lemma 2 Assume that 7, £ F. Then 

\D(p 7 ,p Y )\ < (2h)- 1 \\ 1 \\ 1 [|||7 , |||f + (a* + B»)||7 , || 

1^(7,701 <D(p hl ,p lYl ) . 

Proof. Expanding 7 and 7' in eigcnfunctions, see Eq. ([6]), we get by the Schwarz 
inequality 



(10) 
(11) 



dx dy 



E 



dxdy 



< 



EJA M IM*)l 2 £JA'Jie(y)l* 



|x-y 

Pi 7 i( x )pi7'i(y) 



dxdy 



x- 



dx dy 



which shows Eq. (jlip . To prove Eq. (|10[) . we remark that by Hilbcrt's inequality 
in two dimensions, i.e., |V| > ft.|a;| _1 with h := 47r 2 /(r(l/4)) 4 , and the diamagnetic 
inequality [221, we have 



<M£ M (aO| 5 



iV + A |O. 



(12) 



The claimed bound follows now by multiplication with | A M | and summation over fi 
and v. □ 
We note that because of Eq. (|10p and Lemma [TJ 

Lemma 3 For 7 £ F, we have 

(2hy 1 \\ 1 \\ 1 [iHI*. + ( At + i3 1/2 )||7||i] > Q[|7|] > • 

In a mean-field picture, relativistic electrons are described by one-particle density 
matrices 7 with certain additional properties. In particular, since electrons are 
fcrmions, they obey the Pauli principle and cannot occupy states in the Dirac sea 
which is given by the negative spectral subspace of a one-particle Dirac operator Dj± 
with electromagnetic vector potential A. Mathematically, this is implemented by 
requiring that < 7 < := X(o,oo)(Aa)- As indicated already above, we will choose 
A := A, a choice known as the Furry picture. It is then useful to introduce several 
sets of one-particle density matrices 7 for the subsequent discussion. 

Definition 2 We define the following sets of one-particle density matrices for given 
(maximal) particle number q £ M + 

:= {7 G F I < 7 < A+} , 

:= {7 G I < tr( 7 ) < q} 

4f : = {7 G I tr( 7 ) = q} . 

We note that all sets are closed subsets of F. 
convex. They are only of technical importance, whereas we are ultimately interested 
in describing a system with a fixed number of electrons q, i.e., in minimizing the 

(A) 

energy over the set Sg q , eventually for the quantized case q = N £ N. In physical 



(13) 
(14) 
(15) 

Furthermore, the first two are 
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terms, the projection 1 — can be interpreted as the one-particle density matrix of 
the Dirac sea which we consider frozen. The energy of such a system in Hartrce-Fock 
approximation is given by the functional £ : F — > R as 

S( 1 )=tr[(D A -nh] + aQ[ 1 ]. (16) 

The corresponding rclativistic model has been successfully used in atomic and 
molecular physics [121 HH H7J QI5] . This projection approach implies that all negative- 
energy states and the zero modes are occupied, and N additional particles are then 
added on top of the filled Dirac sea, see also Ref. [11] . 

We now address the mathematical properties of £ (7) in a magnetically confined 
graphene nanostructure. 

Lemma 4 The energy functional £ is well defined and continuous in the \\ ■ \\p norm. 
Furthermore, £ (7) > — fiq for 7 G Sq . 

Proof. The first two claims follow from the definition of the norm, Lemma [3] The 
lower bound is immediate since the only negative term is —fitvy. □ 

Lemma 5 The energy functional £ is coercive on Sq A \ i.e., £("f n ) 00 ifj n € S q A ^ 
and ||7„||f -> 00. 

Proof. For tp in the domain of Da, which equals the one of p Q — fx, we have because 
of the relative compactness of the perturbing term P that 

\\D A i>\\ > W - a*MI - ll^ll > (! - e )ll^o - - mu\\ 

for an arbitrarily small positive e and some MeR. Thus, squaring the inequality and 
taking operator square roots, we get \Da\ > °i\p — A 1 ! — c 2 for some positive constants 
ci and c 2 . Thus for 7 € S( A \ 

tr[(D A - /i) 7 ] = tr((|£> A | - ^) 7 ) > Cl \\>y\\ F - (c 2 + (i)q , 

which implies the coercivity, since Q is non- negative on S^ A \ □ 
In order to fulfill the relative compactness requirement Q, it is in fact sufficient 
to show relative compactness with respect to the free Weyl operator (without the 
magnetic field). Since the following discussion does not rely on this result, it has been 
relegated to the Appendix. 



3.2. Minimization of the Energy 

We now follow Barbaroux et al. |27| and wish to show the existence of a minimizer for 
the Hartree-Fock energy functional (|16|) . Here we will consider only particle numbers 
q which are so small that always 

MS(SW) =in£e(s$) , (17) 

i.e., the magnetic dot is not yet saturated with electrons. Our proof strategy is 
standard: First we show that it is enough to minimize over density matrices of 
finite rank. Then we show that there is a projection with the same particle number 
yielding a lower energy unless it is itself a projection. Moreover, the particle number 
is automatically quantized so that it is enough to require ggN. Equation (|17p thus 
corresponds to the case N < N c discussed in the Introduction. Finally, the existence 
of a minimizer follows by a compactness argument. Let us now go through the steps 
of the proof. 



Multiparticle equations for graphene 



8 



3.2.1. Reduction to Density Matrices with Finite Spectrum. 

Lemma 6 Assume 7 G ^dg • Then there exists a sequence of finite rank density 
matrices Tr- G S 9q such that \\"/k — 7||f — > as A" — > 00. 

Proof. Let fc £ N be a complete set of eigenfunctions of 7. If all eigenvalues are or 
1, the claim is immediate since then 7 is of finite rank itself, and 7 is trace class. Thus 
we can assume that there is an eigenvalue A„ G (0, 1). Now set ck '■= q — SfcLi ^fc- 
Then ck is a non-negative monotone decreasing sequence tending to zero. Define 
Ik ■= X)k=i^fcl£fc)(£fc| + eK\£n){£n\- We now assume n < K and K so big that 

(A) 

A n + £k < 1. Obviously, < Jk € Sa« an( ^ eacn 7^ nas nn ite rank. We now show 
that 7^ — >■ 7 in the _F-norm as A" -> oo. We have 

oo 

fc=JST+l 

Thus we obtain 

oo 

\h-7K\\F< ^tT(\p-fl\\^ k )^ k \) + e K (U,\t>-^n) ■ 

k=K+l 

The first term tends to zero since |ri — /i| 1//2 7|ji — n] 1 ^ 2 G &i(F)a), and the second 
tends to zero since ck — > 0. □ 
The following is an immediate consequence of the continuity of £ in the A-norm 
and the preceding density result: 

Corollary 1 Assume that q > 0. Then 

MS, (4f ) = inf {^( 7 )|7 G S$, rank( 7 ) < ^} . 

3.2.2. Reduction to Projection. In the following we assume q > 0. As customary, 
[q] := max{fc G Z|fc < q} denotes the integer part of g, and we set e q := q — [q]. 
Following the lines of Bach [25], we get: 

Lemma 7 Assume < 7 G Sg^ with finite rank. Then there exists a projection 
A G Sg[q] o,nd a self- adjoint rank one operator R with AR = and tr R = e q such that 

S(A + R) <£(j) . 
Equality holds only if 7 is already of that form. 

Proof. Suppose that 7 is not of that form. Then there exist at least two eigenvalues 
A, A' G (0, 1) of 7; we denote the corresponding normalized eigenvectors by u and 
v. We set 7 := 7 + eS, where S := \u)(u\ — \v)(v\. Note that 7 G Sg q as long as 
< A + e < 1 and < A' — e < 1, which is the case for e in a neighborhood of zero. 
We get 

£(7) - £(7) = e[tr(D A S) + 25RQ( 7 , S)} + e 2 Q(S, S) . 
By explicit computation and use of the Schwarz inequality, we find Q(S, S) < 0, since 
S is a difference of two orthogonal rank one projections. Now - depending on the sign 
of the coefficient linear in e - we lower the energy by increasing or decreasing e from 
zero, until one of the constraints < A + e, A' — e < 1 forbids any further increase 
or decrease of e. This process leaves all the eigenvalues of 7 unchanged except for 
A and A', one of which becomes either or 1. Since there are only finitely many 
eigenvalues of 7 strictly between zero and one, even if they are counted according 
to their multiplicity, iterating this process eliminates all eigenvalues that are strictly 
between and 1, i.e., we have found a density matrix A such that A 2 = A. □ 
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3.2.3. Criterion for Maximal Charge 

Lemma 8 Assume that for 7 <E Sq A ^ with tvy < q, the operator A\h^A^ has at 
least q negative eigenvalues. Then 

If in addition < 7 is a minimizer of £ in S q A \ it follows that try = q. 

Proof. That the left side bounds the right side from above is obvious. To prove the 
converse inequality, we assume that < 7 £ Sq with try < q. By Lemma [3 we can 
assume that 7 is a projection A plus a rank one operator. In particular, the dimension 
of the range of A is at most [tr 7] . Since the dimension of the discrete spectral subspacc 
X of A^/i^pA^ is larger than q, we can find u £ X(l A(SjA) ± with ||u|| < 1 and define 
7 := 7 + u; with u> := \u)(u\. We then get 

5( 7 + u) - £ (7) = \x\(D A - n)u] + 23tQ( 7 , u) = (u, h { ^u) < . 

Therefore this construction yields a density matrix 7 with strictly smaller energy and 
a trace that can be made larger by min{l, q — try}. Iteration of the construction yields 
the desired result. This proves both claims. □ 



3.2.4- Existence of a Minimizer We now wish to show the existence of a minimizer 
by weak lower semi-continuity of the functional on a minimizing sequence and weak 
compactness. This has been addressed by Licb and Simon in the context of orbitals 
in the non-relativistic setting. For density matrices, it was addressed by Solovej [3D] 
in the non-relativistic context, and by Barbaroux et al. [37] for relativistic systems. 

Theorem 1 Assume < q and let \x be in the intersection of the resolvent set of Da 
and the interval (0, where l\ is the first positive eigenvalue ofp Q . Then there exists 

a 7 £ S q A ^ such that 

£( 7 ) = mf£(sW) . 
Moreover, 7 = A + with A a projection, A£ = 0, and ||£|| < 1. 

Proof. Let 7 „ be a minimizing sequence in S q A \ i.e., £("f n ) converges to inf £ (j5 q A ^ . 

Because of the coercivity of £ on Sq A ^ (Lemma [SJ, the sequence 7 „ is bounded in F. 
Thus, according to Banach and Alaoglu, 7n - if necessary by extracting a subsequence 
- converges in the weak-* topology, i.e., there exist 700 £ F such that for all compact 
K , we have 

tr(A^ - M| 1/2 7n|? - Ml 172 ) tv(K\p Q - /*| 1/2 7oo|? - m| 1/2 )- 

Since 

(V>,7ooV>) = 

= tr(|j* - ri-WWmto V\- 1/2 \P H 1/2 7oo||* - mi 1/2 ) 

= Jdm trd^ - mI-V^X^I^ - ^ 2 ln \$ - m| 1/2 ) 

= lim (tp,y n ip) > , 
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we have > 0. Similarly, 7^ < Aj[. Picking an orthonormal basis ei, e 2 , . . ., Fatou's 
lemma gives - possibly under extraction of yet another subsequence - 

q > lim tT7„ > liminf tr(|e„) {e„|7 n ) = . 

n— >oo * — ■* n— >oo 

v 

(A) 

Thus the trace may only decrease. Joining these results shows that Sq is weakly-* 
closed, i.e., 700 € S^. We now show lower semi-continuity in the weak-* topology. 
Concerning the one-particle part, we set Aq" := X(-oo.o](^ ) an( i compute, for some 
i? > and smaller than the first positive eigenvalue of Da, 



/ ^ A o(? -V + P + iV^Ptyo ~ * + if?)" 1 , (18) 



where we used [3TJ Lemma 5.6] sgniJ = tt 1 {irj + H) 1 drj, meant as the Cauchy 
principal value 

r+r 



1 /' +r 

sgn(iJ) = — lim / (i?y + iJ) -1 d7? 

7T r-^oo y 



in the strong topology. Equation (|18j) shows that the product of these two projections 
is compact, expressing that the orthogonality of the positive and negative spectral 
subspaces is not perturbed too much by P, see Eq. Qj. Moreover, it is easy to see 
that Ka '■= \p — A t | 1 ^ 2 A(7Aj l |ji — [i\~ x l 2 is also compact. We are now in a position 
to show the lower semi-continuity of the one-particle part: 

liminftr[( J D A -/i)7„] = (19) 

n—too 

liminf tr ty Q - /i) 7 „ + \p Q - ii\~$P\p - mI~% ~ H\hn\p - lA 
= Uminftr(| | * - M | 1 /V| F k -M| 1/2N 

- 2tr (k a Ka\% - ^r /2 7oo|^ - M| 1/2 ) + tr(P 7 oo) . 

Now, using Fatou's lemma and picking an orthonormal basis ei,e2, . . ., we have for 
the first summand 

hrninftr(|^ - M |VVij?i -M| 1/2 ) 

= [l e -X e -II^O ~ V\ 1/2 ln\$ ~ V\ 1/2 

V 

^ J2 l[ ™™ itT [l e -)( e -l^o ~lA 1/2 ln\p -^| 1/2 ] = booll-F • 

v 

Inserting this into the last line of Eq. (|19[) and undoing the first steps again gives the 
desired bound liminf„_j. 0O tr[(D_A — fJ.)j n ] > tr[(D_4 — /z)7oo]- H remains to show the 
lower semi- continuity of the interaction part which is quadratic in the density matrix. 
Although it is quadratic and positive on S q A \ it is not a positive quadratic form on 
a vector space, i.e., its lower semi-continuity does not follow immediately from the 
Schwarz inequality. However, we can proceed as follows. First, we note that there is 
some constant C such that C > ||7n||F = J2 v (£v ™ > \p — where the arc 

eigenfunctions of y„ with || < 1. Because of Corallary [T] and Lemma [7J we can 
assume that the sum contains at most [q] + 1 summands. Thus, by Lemma [TU1 also the 
standard Sobolev norm of X& is bounded uniformly in n (and v) for any function 
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X- Thus, possibly by extracting another subsequence, we can assume that the £u X 
converge weakly in the i/ 1 / 2 -norm to £1, . . . , Thus, $ (%) — > £n( x ) almost 

everywhere pointwise, see, e.g., Rcf. [321 Theorem 16.1]. The pointwisc convergence 
allows us to use Fatou's lemma to show that 

liminf Q[ 7n ] > g[|6)(ai + ■ • • + iwxwi] • 

We note that the ways of taking the limits - pointwise and in the weak-* sense - agree, 
i.e., the pointwise limit of ^ v ^j l \x)£}j , '\y) is an integral kernel of [27]. □ 
The critical particle number iV c for which the minimizer 7 on has trace N is 
certainly positive. Unfortunately, for an arbitrary perturbation P, we have not been 
able to obtain effective bounds neither from below nor from above for N c . In the next 
section, we therefore consider a specific choice for the perturbation P, first proposed 
in Ref. [TU], and discuss our results for the numerical solution of the self-consistent 
Hartree-Fock equations. 



4. Numerical results 



In this section, we consider a circularly symmetric magnetic dot in graphene, defined 
by a homogeneous background field B with vector potential Ao(x) = (B/2)(— X2, %i) 
and a purely magnetic perturbation P = <j ■ Ai with 

Ai(x) = -{Br/2)[Q{R r) + {R/r) 2 Q{r R)] ( \ , (20) 

where x\ = rcos</>, x-i = rsin0, and O is the Heaviside function. This choice implies 
that the magnetic field vanishes inside a disc of radius R around the origin, while 
outside this disc the field is constant and given by B. On a semiclassical level, one can 
expect bound states inside the disc, built by combining the plane- wave solutions inside 
the disc with the cyclotron orbit solutions outside the disc. The existence of bound 
states is shown below under a fully quantum mechanical approach. Note that single- 
particle states centered far away from this disc ("magnetic dot") represent spatially 
localized Landau cyclotron orbits. The dimensionless "missing flux" S := R 2 /2l 2 3 , 
where Ib ■= (c/leB]) 1 / 2 is the magnetic length, then controls the number of bound 
positive-energy single-particle states energetically below the first bulk Landau level. In 
what follows, all energies are given in units of the first Landau level energy, B ■ 

We solve the interacting problem (a > 0) numerically within the self-consistent 
Hartree-Fock approximation discussed in Sec.[3l with density matrices 7 with tr( 7 ) = 
TV constrained to the Hilbert space spanned by positive single-particle energies. The 
familiar 5 = Landau level energies [2], e s := <J\Jn + (j + l/2)8(j), are expressed in 
terms of the set s := (j, n, a) of quantum numbers j 6 Z + 1/2 (angular momentum), 
11 G No (radial index), and the conduction/ valence band index a = ±. (For n = 
and j < 0, only a = — is allowed and spans the zero-energy level e = 0.) The 
corresponding single-particle cigenspinor |s) has the spatial representation [lOj 

*„ „„ ) (f,«:=(x|,) = ^(. mW 4^ ) ), (21) 

where £ := r 2 /2l 2 B is a dimensionless radial coordinate, and we have (s\s') = S ss > 
with J °° d£[(ipnj) 2 + (VQ) 2 ] = 1- Using the generalized Laguerre polynomials L^, the 



Mul 




Figure 1. Single-particle spectrum of the circular magnetic dot l|20[l vs 
dimcnsionlcss missing flux 6. The lowest five positive-energy solutions E a are 
shown; E a is given in units of the first Landau level energy. 

Landau states (|2"Tj) contain the real-valued functions 

= A% e"«/ 2 L^H^O , (22) 

with normalization factors 

_ / (n-e(-j))! _ / n! 

j " V 2 ( n - e (-^ + l^5l) ! ' J " ™V2(n+|j + i|)r 

For n = and j < 0, we have Aj n = and A~ n has to be multiplied by y/2. 

It is then straightforward to express the perturbation P in this basis, and to 
diagonalize the full single-particle problem numerically. Note that j is still a good 
quantum number, and we label the positive energy solutions \a) with energy E a by 
a = (j, k) with k £ N. Numerical diagonalization of the single-particle Hamiltonian 
yields the orthogonal matrix A in the expansion \a) = A s _ a \s). Assuming an inert 
filled Dirac sea, we only keep states with E a > in what follows. The zero-energy 
states are thus included in the filled Dirac sea, i.e., the chemical potential is assumed 
to be just above zero. The energies E a are shown as a function of the missing flux 
5 in Fig. [TJ Bound states correspond to states with energy < E a < 1 that are 
localized near the origin. As shown in Rcf. [321, t nere are only finitely many states 
below /i - note that /x as defined after Eq. ([SJ) is a constant in (0,1) which we have 
not yet specified - whereas there are infinitely many states with energy between \i and 
1. These states with energy close to 1 are not localized near the dot and correspond 
to weakly perturbed Landau states far away from the dot. We will now pick fj, such 
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Figure 2. Critical particle number N c vs interaction strength a for three different 
values of the missing flux <5. The symbols represent the computed values, and lines 
connecting them are guides to the eye only. 



that it implements this intuition. We have checked that under the choice /i = 0.99, 
all states with < E a < \i correspond to bound states localized near the dot while 
those with fi < E a < 1 correspond to states far away from the dot. From Fig. [TJ we 
can then read off the number of bound states within the magnetic dot. 

Next we address the interacting multiparticle problem, where N electrons are 
added on top of the filled Dirac sea. With the numerically obtained matrix A, the 
two-particle interaction matrix elements 

Va 1 a 2 a 3 a i '■= ^ A Sl )(Jl A S2 t a 2 V Sl S2 s 3 Si A S3 : a 3 A Sl >a4 

Sl,S2,S3,S4 

follow from the Landau-state matrix elements 

T/ S1S2S3S4 =^=j jf^7| (*I 4 • * 8l )(x) • *. a )(x') , (23) 

where s, = (ji,rii,<Ji) and lengths (energies) are expressed in units of Ib (the first 
Landau level energy). Angular momentum conservation dictates ji -\-j2 = js + J4, and 
numerically we find that only momentum exchange processes with k := | j'4 — ji\ < 4 
need to be kept, cf. also Ref. [11]. A useful but lengthy explicit expression for the 
matrix elements (|23j) is given in |Appcndix B[ 

The Hartree-Fock scheme to determine the ground state energy Em for N particles 
with given interaction strength a and missing flux 6 is then standard [HI [22] . 
Numerical calculations were carried out by restricting the Hilbert space to — 18 < j < 2 
and n = 0,1,2, which spans the relevant low-energy sector and very accurately 
describes all stable (i.e., N < N c ) multiparticle ground-state energies reported below. 

The self-consistent numerical solution for the density matrix 7 also allows to 
read off the number of bound electrons in the interacting dot. In particular, only a 
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Figure 3. Hartrcc-Fock ground state energy En vs N for several a with <5 = 3. 
Inset: Same for 5 = 1. 



maximum number N c — N c (a, S) of electrons can be bound by the magnetic dot. and 
for N > N c , we find that N — N c electrons enter bulk Landau states centered far away 
from the dot (and from each other) in order to minimize the Coulomb interaction 
energy. The diagonal elements 7 aa yield the occupation probability of state \a), and 
we can thereby directly infer N c (a, S) from the converged density matrix. The result 
is shown in Fig. [2] for several values of S. Note that N c (6) for a = follows directly 
from Fig.rfJ When increasing the repulsive interaction strength a, electrons tend to be 
pushed out of the dot, and N c monotonically decreases with a. For sufficiently strong 
interactions, a > 1.75, and moderate values of the missing flux, 6 < 3, we find that 
only a single electron can be bound by such a dot (N c = 1). Coulomb interactions in 
graphene nanostructures are thus very significant for the physically relevant regime 
a < 2. 

Hartree-Fock results for the ground state energy En are shown for several S and 
a in Fig. [31 For N > N c , the shown energies still depend on the chosen Hilbert space 
dimension, and for bigger basis size, they move to smaller values. This is clear on 
physical grounds, since N — N c electrons will stay as far away as possible from the dot 
region and from each other in order to minimize the Coulomb energy. In the limit of 
infinite basis size, the Coulomb interaction energy can be minimized by moving N — N c 
electrons to infinity, such that E^ + i — En = 1 for N > N c . Obviously, for the chosen 
basis size, the results in Fig. [3] do not yet match onto this asymptotic form. However, 
the shown results for Em with N < N c in Fig. [3l describing the stable multiparticlc 
case, are accurate and do not change when increasing the basis size. This behavior 
gives an additional criterion to find the number iV c and thereby the sought stability 
condition. 
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5. Discussion 

In this paper, we have studied the interacting multiparticle problem in graphene 
quantum dots. Such quantum dots can be created in a tunable way by imposing 
suitable inhomogeneous magnetic field profiles. After providing the mathematical 
foundations of Hartree-Fock theory and discussing the conditions for the existence of a 
minimizer, we have given predictions for the maximum number of bound electrons in a 
specific example. The Hartree-Fock ground state energies shown here represent upper 
bounds for the true ground state energy. An alternative to Hartree-Fock calculations 
is to employ the Miiller functional [Ml Hi]. As opposed to the Hartree-Fock functional, 
the exchange energy is no longer dominated by -D[p 7 ] but instead by the kinetic energy. 
In fact, the Miiller correlation energy shares the feature of the Dirac exchange energy 
ptj z that it underestimates the correlation energy, resulting in too low energies. We 
expect that the Miiller functional then implies lower bounds for the ground-state 
energy, and we can thereby get both upper and lower bounds for the exact result. 
This work is currently in progress. 
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Appendix A. Some useful lemmata 

The relevance of the following result is that it suffices to show relative (form) 
compactness with respect to the free Weyl operator to fulfill the compactness 
requirements (f4]). 

For some R > 1 (not to be confused with the radius R used in Sec. 0}, we define 
a smooth cutoff function xr & C^°(R 2 ; [0,1]) with the property that \ R = 1 for 
|x| < R/2 and X R = for |x| > R; it follows that ||Vxjj||oo = 0(1/ R). We will often 
use that for tp e (M 2 : C 2 ) , 

[Da,Xr]<P = -Vxh<^ • (A.l) 

Lemma 9 Let A = (A\, A%) be a magnetic vector potential with Aj £ Lf^ c (R 2 ) for 
j = 1,2. Let T be a bounded symmetric (matrix-valued) multiplication operator such 
that 

a) ||T(x)|| -> as ||x|| -> oo, 

b) T\tr ■ p + i| -1 / 2 is a compact operator. 

Then, for any A in the resolvent set of Da, the operator T\Da + A| -1 / 2 is compact. 

Proof. Clearly the claim follows if we prove that xrT\D a + i|~ 1//2 is compact for all 
R > 1, since xrT\Da + i|~ 1//2 — s- T\D a + i| -1 ^ 2 in the operator norm as R — >• oo and 
the operator \Da + ip^-DA + A| -1 / 2 is bounded. We observe that 

XrT\D a + i\- 1/2 - XrTW ■ p + r 1 / 2 + 

+ T\tr ■ p + i|- 1/2 |<T ■ p + i\ 1/2 X R(\D A + r 1/2 ~ W ■ P + ir V2 ) • 
Therefore, it suffices to show that |er • p + i^^XrQDa + i| -1 ^ 2 — \tr ■ p + i|~ 1//2 ) is a 
bounded operator. In order to do so, we first note that the following resolvent identity 
holds in L 2 (K 2 : C 2 ): For z € iK \ {0}, write 

R A (z) = (D A + z)-\ Ro(z) = (cr-p + z)" 1 , 
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then 

X r(Ra(z) - Ro(z)) (A.2) 
= Ro(z)(ia ■ W X r)(Ra(z) - R (z)) - R (z)a ■ A X rRa(z) 
=: R (z)T(z) . 

Before proving the above equation, using the spectral theorem, we compute 

dt 1 

1 r 00 at 



\D A +i\-v*=(Di+i)-v*= 1 r 



[^(-ivT+t) - R A {iVT+t)] 



V2tt Jo 2it l /4 1 /TTt 

An analogous formula holds for |cr ■ p + i| -1 / 2 with R A replaced by Rq. Therefore, 
using also Eq. (|A~2l) we get 

w ■ P + (\d a + r 1 ' 2 k • p + r 1/2 ' 
dt 



- v kI 



|cr p + i| 1/2 i? («iVl + t)T(K,iy/l + t) 



Noting that \tr ■ p + i\ l/2 R a {ni^/TTt) is bounded and that ||r(«;i\/l + i)|| < c/y/l + t 
for some constant c, we conclude that the operator above is bounded. 

It remains to show Eq. (|A~2|) . For 4> € C£°(R 2 : C 2 ), using Eq. (jA~T) we find 

{R A (z) - R {z))xr{(t ■ p + z)4> 

= (R A {z) - Ro{z))[icr ■ V X r + (<r ■ P + z)xr]</> 

= (i? j4 (z) - i?o(^))i<T • Vxh^ - i?A(^)cr • Axr4> 

= [{Ra{z) - Ro(z))icr ■ ^xrRo(z) 

- R A (z)(T ■ AxrRo{z)]((t -p + z)<f>. 

Since the range of (cr ■ p + z)|Cg°(R 2 : C 2 ) is dense in L 2 (R 2 : C 2 ), we obtain the 
adjoint of Eq. (|A.2|) by a limiting argument. □ 
The next result allows us to use weak H 1 / 2 convergence in the proof of the main 
theorem. 

Lemma 10 Assume for the magnetic field and the magnetic vector potential that their 
components are locally bounded, i.e., Aj,Bj G ij^ c (R 2 ) for j = 1,2. Then, for any 
X G Cg°(R 2 ; [0, 1]), there exists a constant c x > such that for all <p G C£°(R 2 : C 2 ), 
we have 

(0, X|P|X0) < 2 3/2 (0, \? - /#) + c x ||0|| 2 . (A.3) 

Proof. A simple application of the inequality (a + b) 2 > (1 — <5)a 2 + (1 — l/<5)6 2 yields, 
as quadratic forms on C^(R 2 : C 2 ), 

- v) 2 x > (1 - %[(p + a) 2 + cr • b] x + (1 - r Vx a 

> (i-i)Vx+(i-*)(i-OxA J x 
+ (i-5)x Cr -B X + (i-r 1 ) M V 

> \xV 2 X - c 2 , 
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where in the last inequality we set S — 1/2 and c 2 = ||xA 2 x||oc + ||x|B|x||oo/2 + /i 2 . 
With the help of the above inequality and Eq. (|A. 1|) , we obtain, for any <j> € C§° (R 2 : 

C 2 ), 

\{<t>> \x\p\x? d>) < i||| P |x0ll 2 < ||(? - m)x0II 2 + C 2 II^II 2 

<2||(^- M )^|| 2 + 2|| ( t.Vx^|| 2 + c 2 ||0|| 2 

< 2(0, (? - M ) 2 0) + (0, [2||Vx||oo + c 2 ]0) . 

Using this last estimate and the fact that the square root is operator monotone, we 
get Eq. (pO} with c x = 2[2||Vx||oo + c 2 ] 1 / 2 . □ 

Appendix B. Two-particle matrix elements 

In this appendix, we provide the explicit form of the two-particle Coulomb matrix 
elements ([231) in terms of the Landau level spinors (|2"Tj) . With the radial functions 
(|22[) and momentum exchange k = \j& — we have 

/■7T POO /"O0 /] i\ 

a / / / , cos(k<p) 



where V'l ■ •= ^Pt^Pt + &ic r 2'4'i' t p2 ■ We now employ the expansion formula 

1 1 00 

7— =4^ 2 = ^=E(W£>) £/2 ^(cos<« 



with Legendre functions P^(cos0), where £> = max(^(') and = min(£,£'). The 
angular integration can then be performed by using the relation 

^ cos(fc<ft)Pi(cos0) = 2 ^ 2 ip " C(e + k)/2- t t ■ (B.2) 



7T 



The coefficients C m; £ with m G No are the expansion coefficients of a hypergeometric 
function, and C(£ +k y 2 -,e om y f° r even k + I and £ > k. In particular, Co : £ = 1, 
while for < m < £, we have the product representation 

_fV (z-l/2)(£+l-Q 
W-ll i(f + 1 / 2 _j) • 

i— 1 x ' ' 

To perform the £, £' integrations in Eq. (jB.ip . we insert the explicit form of in 
Eq. (f2"2"j) . with the generalized Laguerre polynomials (n,m G No) 



z! \ n — i 

After some algebra, with the I summation only extending over I + k G 2Z, we find the 
lengthy result 

V Sl S 2 S 3Si = a E E " 2^ + ^|" C ^+ k )/ 2 ^ ( B ' 3 ) 



x (*„,+ + <y, t _tr l0 - 4 ) (<V,+ + S, r ,-a 2 a 3 ) A^A^A^A* 
H 



Si 

ni n 2 n 3 n 4 / Nmi+rra2+m3+ „ i4 



E E E E 



mi!m 2 !m 3 !m4! 

mi- U 7712 —U 771-3=0 7714— U 
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x 



X 



ni + lii - VM \f n 2 + |ja - r?7 2 l 
ni — mi y y n 2 - m 2 

% + |J3 - V/2| ]^4+ |J4 - »7/2| 

n 3 -m 3 J \ h 4 -m 4 
x J(M + M', Af') + (si, s 4 ) <-> (s 2 , s 3 ) , 
where h 1A := 7i M - 5 J)j+ 6(-ji >4 ), n 2 ,3 := n 2 , 3 - <5 t ,/,+ e (~:/ 2 ,3), and 

M := TO! + 77 l4 + ^ (|ii - f?/2| + |j 4 - ?//2| - , 

M' := to 2 + to 3 + l - (\ J2 - ?//2| + |j 3 - r)'/2\ + t) ■ 
Finally, for 71,77,' € No, we have 

(277 + 1)!! f 1 y n ' 
I{n,n):= / dy 



(1 + 7,)"+ 3 / 2 ' 

This allows for the numerical evaluation of the Coulomb interaction matrix elements, 
since all summations in Eq. (|B.3[) converge rapidly. Finally, note that the matrix 
elements obey the symmetry relations 

Vs 1 s 2 s 3 s i = Vs 2 s 1 s 4 ,s 3 = Ks 4 s 3 s 2 si • (B-4) 
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